#############################################
# differential abundance analysis
#############################################

library(ggordiplots)
library(microViz)
library("vegan")
library(phyloseqCompanion)
library(ggplot2)
library(Rmisc)
library(ggsignif)
library(ggpubr)
library(dplyr)
library("gridExtra")
library(microbiomeMarker)
library(tidyr)
library(FUNGuildR)
library(SciViews)

# set directory
setwd("~/Research/Mycorrhizae/Fieldwork/Genetics/NGS sequencing/Code")

# set scientific notation off
options(scipen=999)

# load phyloseq variables
load("~/Research/Mycorrhizae/Fieldwork/Genetics/NGS sequencing/Code/phyloseq final variables.RData")

# get otu abundance matrix from phyloseq object
all.otu.abundance <- psmelt(physeq)

# form combined column of taxonomy separated by comma
all.otu.abundance$Combined<-paste(all.otu.abundance$Kingdom, all.otu.abundance$Phylum, all.otu.abundance$Class, all.otu.abundance$Order, all.otu.abundance$Family, all.otu.abundance$Genus, all.otu.abundance$Species,sep=",")

# assign guilds
sample_guilds <- funguild_assign(all.otu.abundance$Combined, db = get_funguild_db())

# turn assignments into dataframe
sample_guilds <- as.data.frame(sample_guilds)

# combine with abundance dataframe
all.otu.abundance$trophic.guild <- sample_guilds$trophicMode
all.otu.abundance$confidence <- sample_guilds$confidenceRanking

# now, manually fill in genera that need more specification for trophic mode

# Entoloma (run in this order)
all.otu.abundance$trophic.guild[which(all.otu.abundance$Genus=="Entoloma")] <- c("Saprotroph")
all.otu.abundance$trophic.guild[which(all.otu.abundance$Genus=="Entoloma" & all.otu.abundance$Species=="byssisedum")] <- c("Saprotroph")
all.otu.abundance$trophic.guild[which(all.otu.abundance$Genus=="Entoloma" & all.otu.abundance$Species=="caccabus")] <- c("Saprotroph")
all.otu.abundance$trophic.guild[which(all.otu.abundance$Genus=="Entoloma" & all.otu.abundance$Species=="calobrunneum")] <- c("Saprotroph")
all.otu.abundance$trophic.guild[which(all.otu.abundance$Genus=="Entoloma" & all.otu.abundance$Species=="chamaemori")] <- c("Saprotroph")
all.otu.abundance$trophic.guild[which(all.otu.abundance$Genus=="Entoloma" & all.otu.abundance$Species=="cuspidiferum")] <- c("Saprotroph")
all.otu.abundance$trophic.guild[which(all.otu.abundance$Genus=="Entoloma" & all.otu.abundance$Species=="fernandae")] <- c("Saprotroph")
all.otu.abundance$trophic.guild[which(all.otu.abundance$Genus=="Entoloma" & all.otu.abundance$Species=="formosum")] <- c("Saprotroph")
all.otu.abundance$trophic.guild[which(all.otu.abundance$Genus=="Entoloma" & all.otu.abundance$Species=="fuscotomentosum")] <- c("Saprotroph")
all.otu.abundance$trophic.guild[which(all.otu.abundance$Genus=="Entoloma" & all.otu.abundance$Species=="hirtipes")] <- c("Saprotroph")
all.otu.abundance$trophic.guild[which(all.otu.abundance$Genus=="Entoloma" & all.otu.abundance$Species=="huijsmanii")] <- c("Saprotroph")
all.otu.abundance$trophic.guild[which(all.otu.abundance$Genus=="Entoloma" & all.otu.abundance$Species=="jubatum")] <- c("Saprotroph")
all.otu.abundance$trophic.guild[which(all.otu.abundance$Genus=="Entoloma" & all.otu.abundance$Species=="kristiansenii")] <- c("Saprotroph")
all.otu.abundance$trophic.guild[which(all.otu.abundance$Genus=="Entoloma" & all.otu.abundance$Species=="llimonae")] <- c("Saprotroph")
all.otu.abundance$trophic.guild[which(all.otu.abundance$Genus=="Entoloma" & all.otu.abundance$Species=="luteofuscum")] <- c("Saprotroph")
all.otu.abundance$trophic.guild[which(all.otu.abundance$Genus=="Entoloma" & all.otu.abundance$Species=="mirum")] <- c("Saprotroph")
all.otu.abundance$trophic.guild[which(all.otu.abundance$Genus=="Entoloma" & all.otu.abundance$Species=="nitens")] <- c("Saprotroph")
all.otu.abundance$trophic.guild[which(all.otu.abundance$Genus=="Entoloma" & all.otu.abundance$Species=="palustre")] <- c("Saprotroph")
all.otu.abundance$trophic.guild[which(all.otu.abundance$Genus=="Entoloma" & all.otu.abundance$Species=="proterum")] <- c("Saprotroph")
all.otu.abundance$trophic.guild[which(all.otu.abundance$Genus=="Entoloma" & all.otu.abundance$Species=="rhombisporum")] <- c("Saprotroph")
all.otu.abundance$trophic.guild[which(all.otu.abundance$Genus=="Entoloma" & all.otu.abundance$Species=="sphagneti")] <- c("Saprotroph")
all.otu.abundance$trophic.guild[which(all.otu.abundance$Genus=="Entoloma" & all.otu.abundance$Species=="terreum")] <- c("Saprotroph")
all.otu.abundance$trophic.guild[which(all.otu.abundance$Genus=="Entoloma" & all.otu.abundance$Species=="tjallingiorum")] <- c("Saprotroph")
all.otu.abundance$trophic.guild[which(all.otu.abundance$Genus=="Entoloma" & all.otu.abundance$Species=="vernum")] <- c("Saprotroph")
all.otu.abundance$trophic.guild[which(all.otu.abundance$Genus=="Entoloma" & all.otu.abundance$Species=="zuccherellii")] <- c("Saprotroph")
all.otu.abundance$trophic.guild[which(all.otu.abundance$Genus=="Entoloma" & all.otu.abundance$Species=="sericeum")] <- c("Saprotroph")

all.otu.abundance$trophic.guild[which(all.otu.abundance$Genus=="Entoloma" & all.otu.abundance$Species=="sericatum")] <- c("Symbiotroph")
all.otu.abundance$trophic.guild[which(all.otu.abundance$Genus=="Entoloma" & all.otu.abundance$Species=="borgenii")] <- c("Symbiotroph")
all.otu.abundance$trophic.guild[which(all.otu.abundance$Genus=="Entoloma" & all.otu.abundance$Species=="paludicola")] <- c("Symbiotroph")
all.otu.abundance$trophic.guild[which(all.otu.abundance$Genus=="Entoloma" & all.otu.abundance$Species=="politum")] <- c("Symbiotroph")
all.otu.abundance$trophic.guild[which(all.otu.abundance$Genus=="Entoloma" & all.otu.abundance$Species=="rhodopolium")] <- c("Symbiotroph")
all.otu.abundance$trophic.guild[which(all.otu.abundance$Genus=="Entoloma" & all.otu.abundance$Species=="rubrobasis")] <- c("Symbiotroph")
all.otu.abundance$trophic.guild[which(all.otu.abundance$Genus=="Entoloma" & all.otu.abundance$Species=="serpens")] <- c("Symbiotroph")

# for Oidiodendron, only mauis is mycorrhizal and everything else is saprotrophic (run in this order)
all.otu.abundance$trophic.guild[which(all.otu.abundance$Genus=="Oidiodendron")] <- c("Saprotroph")
all.otu.abundance$trophic.guild[which(all.otu.abundance$Genus=="Oidiodendron" & all.otu.abundance$Species=="maius")] <- c("Symbiotroph")

# for Sistotema, only luteoviride and muscicola are mycorrhizal and everything else is saprotrophic (run in this order)
all.otu.abundance$trophic.guild[which(all.otu.abundance$Genus=="Sistotrema")] <- c("Saprotroph")
all.otu.abundance$trophic.guild[which(all.otu.abundance$Genus=="Sistotrema" & all.otu.abundance$Species=="luteoviride")] <- c("Symbiotroph")
all.otu.abundance$trophic.guild[which(all.otu.abundance$Genus=="Sistotrema" & all.otu.abundance$Species=="muscicola")] <- c("Symbiotroph")

# for Cadophora, orchidicola is ericoid mycorrhizal 
all.otu.abundance$trophic.guild[which(all.otu.abundance$Genus=="Cadophora")] <- c("Saprotroph")
all.otu.abundance$trophic.guild[which(all.otu.abundance$Genus=="Cadophora" & all.otu.abundance$Species=="orchidicola")] <- c("Symbiotroph")

# wasn't sure about Trichophaea, Meliniomyces (some can be but no species information). For genera with unknown species, if saprotroph was primary
# lifestyle in fungal functional traits database, then used that

# label symbiotrophs for those generas that not automatically assigned proper trophic guild from funguild
all.otu.abundance$trophic.guild[which(all.otu.abundance$Genus=="Tomentella" | all.otu.abundance$Genus=="Hygrophorus" | 
   all.otu.abundance$Genus=="Pseudotomentella" | all.otu.abundance$Genus=="Trichophaea" | all.otu.abundance$Genus=="Clavulina" | 
     all.otu.abundance$Genus=="Wilcoxina" | all.otu.abundance$Genus=="Genabea" | all.otu.abundance$Genus=="Rhodoscypha" | 
     all.otu.abundance$Genus=="Pulvinula" | all.otu.abundance$Genus=="Tomentellopsis" | all.otu.abundance$Genus=="Otidea" | 
     all.otu.abundance$Genus=="Geopora" | all.otu.abundance$Genus=="Naucoria" | all.otu.abundance$Genus=="Amanita" | 
     all.otu.abundance$Genus=="Hymenogaster" | all.otu.abundance$Genus=="Xerocomus" | all.otu.abundance$Genus=="Hebeloma" | 
     all.otu.abundance$Genus=="Alnicola" | all.otu.abundance$Genus=="Lyophyllum" | all.otu.abundance$Genus=="Boletus" | 
     all.otu.abundance$Genus=="Tricholoma")] <- c("Symbiotroph")

# saprotrophs
all.otu.abundance$trophic.guild[which(all.otu.abundance$Genus=="Mortierella" | all.otu.abundance$Genus=="Solicoccozyma" | 
  all.otu.abundance$Genus=="Thelephora" | all.otu.abundance$Genus=="Hygrocybe" | all.otu.abundance$Genus=="Clavaria" | 
    all.otu.abundance$Genus=="Hygrocybe" | all.otu.abundance$Genus=="Serendipita" | all.otu.abundance$Genus=="Meliniomyces" | 
    all.otu.abundance$Genus=="Clavulinopsis" | all.otu.abundance$Genus=="Chalara" | all.otu.abundance$Genus=="Ramariopsis" |
    all.otu.abundance$Genus=="Glarea" | all.otu.abundance$Genus=="Hymenoscyphus" | all.otu.abundance$Genus=="Acephala" | 
    all.otu.abundance$Genus=="Mycosymbioces" | all.otu.abundance$Genus=="Phialocephala" | all.otu.abundance$Genus=="Infundichalara" |
    all.otu.abundance$Genus=="Odontia" | all.otu.abundance$Genus=="Cudoniella" | all.otu.abundance$Genus=="Thelebolus" |
    all.otu.abundance$Genus=="Crocicreas" | all.otu.abundance$Genus=="Cuphophyllus" | all.otu.abundance$Genus=="Cyathicula" | 
    all.otu.abundance$Genus=="Scutellinia" | all.otu.abundance$Genus=="Pseudaleuria" | all.otu.abundance$Genus=="Xenopolyscytalum" | 
    all.otu.abundance$Genus=="Articulospora" | all.otu.abundance$Genus=="Vibrissea" | all.otu.abundance$Genus=="Ascocoryne" | 
    all.otu.abundance$Genus=="Claussenomyces" | all.otu.abundance$Genus=="Curvibasidium" | all.otu.abundance$Genus=="Chromosera" |
    all.otu.abundance$Genus=="Endogone" | all.otu.abundance$Genus=="Xanthoporia" | all.otu.abundance$Genus=="Pachyella" |
    all.otu.abundance$Genus=="Filosporella" | all.otu.abundance$Genus=="Lamprospora" | all.otu.abundance$Genus=="Hyphodontiella" | 
    all.otu.abundance$Genus=="Triposporium" | all.otu.abundance$Genus=="Onnia" | all.otu.abundance$Genus=="Gliophorus" | 
    all.otu.abundance$Genus=="Byssonectria" | all.otu.abundance$Genus=="Therrya" | all.otu.abundance$Genus=="Taeniospora" | 
    all.otu.abundance$Genus=="Trichocladium" | all.otu.abundance$Genus=="Entocybe" | all.otu.abundance$Genus=="Paraconiothyrium" |
    all.otu.abundance$Genus=="Athelia" | all.otu.abundance$Genus=="Typhula" | all.otu.abundance$Genus=="Cystoderma" | 
    all.otu.abundance$Genus=="Lecythophora" | all.otu.abundance$Genus=="Coniochaeta" | all.otu.abundance$Genus=="Sclerococcum" |
    all.otu.abundance$Genus=="Pseudocosmospora" | all.otu.abundance$Genus=="Coccomyces" | all.otu.abundance$Genus=="Ceratorhiza" | 
    all.otu.abundance$Genus=="Fusicolla" | all.otu.abundance$Genus=="Didymosphaeria" | all.otu.abundance$Genus=="Tephrocybe" | 
    all.otu.abundance$Genus=="Ceratobasidium" | all.otu.abundance$Genus=="Basidiobolus" | all.otu.abundance$Genus=="Fusarium" | 
    all.otu.abundance$Genus=="Aphanobasidium" | all.otu.abundance$Genus=="Cladosporium" | all.otu.abundance$Genus=="Hypsizygus" | 
    all.otu.abundance$Genus=="Pezoloma" | all.otu.abundance$Genus=="Sagaranella" | all.otu.abundance$Genus=="Bjerkandera" | 
    all.otu.abundance$Genus=="Colpoma" | all.otu.abundance$Genus=="Clitopilus" | all.otu.abundance$Genus=="Pseudofabraea" | 
    all.otu.abundance$Genus=="Hypholoma" | all.otu.abundance$Genus=="Rhodotorula" | all.otu.abundance$Genus=="Pluteus" | 
    all.otu.abundance$Genus=="Strobilurus" | all.otu.abundance$Genus=="Pholiota" | all.otu.abundance$Genus=="Psilocybe" | 
    all.otu.abundance$Genus=="Exophiala" | all.otu.abundance$Genus=="Galerina" | all.otu.abundance$Genus=="Pyxidiophora" | 
    all.otu.abundance$Genus=="Scopuloides" | all.otu.abundance$Genus=="Eleutheromyces" | all.otu.abundance$Genus=="Ganoderma" | 
    all.otu.abundance$Genus=="Parafabraea" | all.otu.abundance$Genus=="Armillaria" | all.otu.abundance$Genus=="Ceraceomyces" | 
    all.otu.abundance$Genus=="Mollisia" | all.otu.abundance$Genus=="Septonema" | all.otu.abundance$Genus=="Kuehneromyces" | 
    all.otu.abundance$Genus=="Mycosphaerella" | all.otu.abundance$Genus=="Acrodontium" | all.otu.abundance$Genus=="Halokirschsteiniothelia" | 
    all.otu.abundance$Genus=="Tritirachium" | all.otu.abundance$Genus=="Cladophialophora" | all.otu.abundance$Genus=="Tyromyces" | 
    all.otu.abundance$Genus=="Stropharia" | all.otu.abundance$Genus=="Rhodosporidiobolus" | all.otu.abundance$Genus=="Deconica" | 
    all.otu.abundance$Genus=="Trimmatostroma" | all.otu.abundance$Genus=="Nematoctonus" | all.otu.abundance$Genus=="Juncaceicola" | 
    all.otu.abundance$Genus=="Fomitopsis" | all.otu.abundance$Genus=="Tubaria" | all.otu.abundance$Genus=="Pezicula" | 
    all.otu.abundance$Genus=="Capronia" | all.otu.abundance$Genus=="Phaeosphaeria" | all.otu.abundance$Genus=="Neosetophoma" | 
    all.otu.abundance$Genus=="Postia" | all.otu.abundance$Genus=="Atractiella" | all.otu.abundance$Genus=="Biappendiculispora" | 
    all.otu.abundance$Genus=="Rhexocercosporidium" | all.otu.abundance$Genus=="Pyrenochaetopsis" | all.otu.abundance$Genus=="Coniophora" | 
    all.otu.abundance$Genus=="Clitocybe" | all.otu.abundance$Genus=="Arrhenia" | all.otu.abundance$Genus=="Mycena" | all.otu.abundance$Genus=="Lepista" | 
    all.otu.abundance$Genus=="Fayodia" | all.otu.abundance$Genus=="Melanoleuca" | all.otu.abundance$Genus=="Squamanita" | 
    all.otu.abundance$Genus=="Ripartites" | all.otu.abundance$Genus=="Mycenella" | all.otu.abundance$Genus=="Collybia" | 
    all.otu.abundance$Genus=="Cystodermella")] <- c("Saprotroph")
 
# parasites/pathogens
all.otu.abundance$trophic.guild[which(all.otu.abundance$Genus=="Filobasidiella")] <- c("Pathotroph")
all.otu.abundance$trophic.guild[which(all.otu.abundance$Genus=="Neonectria")] <- c("Pathotroph")
all.otu.abundance$trophic.guild[which(all.otu.abundance$Genus=="Ilyonectria")] <- c("Pathotroph")
all.otu.abundance$trophic.guild[which(all.otu.abundance$Genus=="Volutella")] <- c("Pathotroph")
all.otu.abundance$trophic.guild[which(all.otu.abundance$Genus=="Monodictys")] <- c("Pathotroph")
all.otu.abundance$trophic.guild[which(all.otu.abundance$Genus=="Nectriopsis")] <- c("Pathotroph")
all.otu.abundance$trophic.guild[which(all.otu.abundance$Genus=="Ramularia")] <- c("Pathotroph")
all.otu.abundance$trophic.guild[which(all.otu.abundance$Genus=="Tolypocladium")] <- c("Pathotroph")
all.otu.abundance$trophic.guild[which(all.otu.abundance$Genus=="Metapochonia")] <- c("Pathotroph")
all.otu.abundance$trophic.guild[which(all.otu.abundance$Genus=="Microdochium")] <- c("Pathotroph")
all.otu.abundance$trophic.guild[which(all.otu.abundance$Genus=="Metarhizium")] <- c("Pathotroph")
all.otu.abundance$trophic.guild[which(all.otu.abundance$Genus=="Phacidium")] <- c("Pathotroph")

# check to make sure no genera left in combined trophic guild categories
unique(all.otu.abundance$Family[which(all.otu.abundance$trophic.guild=="Saprotroph-Symbiotroph")])
unique(all.otu.abundance$Genus[which(all.otu.abundance$trophic.guild=="Pathotroph-Saprotroph-Symbiotroph")])
unique(all.otu.abundance$Genus[which(all.otu.abundance$trophic.guild=="Pathotroph-Saprotroph")])
unique(all.otu.abundance$Genus[which(all.otu.abundance$trophic.guild=="Pathotroph-Symbiotroph")])

# we can't look at families because many families have multiple trophic guilds within the family

# make a list of all symbiotrophic species 
all.otu.abundance$names <- c(1:840384)
all.otu.abundance$names <- paste(all.otu.abundance$Genus, all.otu.abundance$Species, sep=" ")
symbio.unique.species <- sort(unique(all.otu.abundance$names[which(all.otu.abundance$trophic.guild=="Symbiotroph")]))
symbio.unique.species <- list(symbio.unique.species)
write.table(as.data.frame(symbio.unique.species),file="mylist.csv", quote=F,sep=",",row.names=F)

# put all labeled saprotrophs together
all.otu.abundance$trophic.guild[which(all.otu.abundance$trophic.guild==" Saprotroph")] <- c("Saprotroph")

# now add exploration type designation

# remove NAs in all otu abundance genera because otherwise if statement won't run
all.otu.abundance.NA <- all.otu.abundance[is.na(all.otu.abundance$Genus),]
all.otu.abundance <- all.otu.abundance[!is.na(all.otu.abundance$Genus),]
dim(all.otu.abundance)

# remove NAs in all.otu.abundance trophic guild category
all.otu.abundance.NA.trophic <- all.otu.abundance[is.na(all.otu.abundance$trophic.guild),]
all.otu.abundance <- all.otu.abundance[!is.na(all.otu.abundance$trophic.guild),]
dim(all.otu.abundance)
dim(all.otu.abundance.NA)
dim(all.otu.abundance.NA.trophic)

# create new column for distance type
all.otu.abundance$distance.type <- as.vector(1:535776)

#create similar column of distance type for NA values of all.otu.abundance
all.otu.abundance.NA$distance.type <- as.vector(1:304416)

#create similar column of distance type for NA values of all.otu.abundance
all.otu.abundance.NA.trophic$distance.type <- as.vector(1:192)

# upload funguild database
funguild <- read.csv("Fungal Functional Traits database.csv", header=TRUE)

# fill in exploration type 
for (i in 1:535776){
  temp.genus <- all.otu.abundance$Genus[i]
  if (all.otu.abundance$trophic.guild[i]=="Symbiotroph"){
    temp.distance <- funguild$Ectomycorrhiza_exploration_type_template[which(funguild$GENUS==temp.genus)]
    all.otu.abundance$distance.type[i] <- temp.distance
  }
  else{all.otu.abundance$distance.type=="NA"}
}

# combine dataframes back together
all.otu.abundance.final <- rbind(all.otu.abundance, all.otu.abundance.NA, all.otu.abundance.NA.trophic)
dim(all.otu.abundance.final)

write.csv(all.otu.abundance.final, "all.otu.abundance.final.csv")

#all.otu.abundance <- read.csv("all.otu.abundance.csv", header=TRUE)

##############################################
# RELATIVE ABUNDANCE CALCULATIONS #
##############################################

# find relative abundance of each trophic guild and distance type by sample

# assign unique number to each sample
library(tidyverse)

all.otu.abundance.final <- all.otu.abundance.final %>%
  group_by(Sample) %>%
  mutate(ID = cur_group_id())

all.otu.abundance.final <- as.data.frame(all.otu.abundance.final)

# new dataset for storing RA of distance types
m <- matrix(0, ncol = 3, nrow = 96)
exploration.guild.RA <- data.frame(m)

# EXPLORATION TYPE #

# calculate RA for each distance type
for (i in 1:96){
  temp <- all.otu.abundance.final[which(all.otu.abundance.final$ID==i),]
  abundance.total <- sum(temp$Abundance)
  abundance.short.delicate <- sum(temp$Abundance[which(temp$distance.type=="short-distance_delicate")])
  abundance.short.contact <- sum(temp$Abundance[which(temp$distance.type=="contact")])
  abundance.short.coarse <- sum(temp$Abundance[which(temp$distance.type=="short-distance_coarse")])
  abundance.short <- (abundance.short.delicate + abundance.short.contact + abundance.short.coarse)/abundance.total
  abundance.medium.fringe <- sum(temp$Abundance[which(temp$distance.type=="medium-distance_fringe")])/abundance.total
  abundance.medium.smooth <- sum(temp$Abundance[which(temp$distance.type=="medium-distance_smooth")])/abundance.total
  abundance.long <- sum(temp$Abundance[which(temp$distance.type=="long-distance")])/abundance.total
  exploration.guild.RA$short.distance[i] <- abundance.short
  exploration.guild.RA$medium.distance.fringe[i] <- abundance.medium.fringe
  exploration.guild.RA$medium.distance.smooth[i] <- abundance.medium.smooth
  exploration.guild.RA$long.distance[i] <- abundance.long
  exploration.guild.RA$Stream[i] <- temp$Stream[1]
  exploration.guild.RA$Distance[i] <- temp$Distance[1]
  exploration.guild.RA$Transect[i] <- temp$Transect[1]
  exploration.guild.RA$Side[i] <- temp$Side[1]
}

# TROPHIC GUILD #

# calculate RA for each trophic guild
for (i in 1:96){
  temp <- all.otu.abundance.final[which(all.otu.abundance.final$ID==i),]
  abundance.total <- sum(temp$Abundance)
  # only include probable and highly probable (remove possible values)
  abundance.saprotroph <- sum(temp$Abundance[which(temp$trophic.guild=="Saprotroph" & temp$confidence != "Possible")])/abundance.total
  abundance.symbiotroph <- sum(temp$Abundance[which(temp$trophic.guild=="Symbiotroph" & temp$confidence != "Possible")])/abundance.total
  exploration.guild.RA$saprotroph[i] <- abundance.saprotroph
  exploration.guild.RA$symbiotroph[i] <- abundance.symbiotroph
}

#########################################################################
# RICHNESS CALCULATIONS FOR TROPHIC GUILD AND EXPLORATION TYPE #
#########################################################################

exploration.guild.RA$saprotroph.richness <- c(1:96)
exploration.guild.RA$symbiotroph.richness <- c(1:96)
exploration.guild.RA$long.richness <- c(1:96)
exploration.guild.RA$short.richness <- c(1:96)
exploration.guild.RA$medium.fringe.richness <- c(1:96)
exploration.guild.RA$medium.smooth.richness <- c(1:96)

for (i in 1:96){
  temp <- all.otu.abundance.final[which(all.otu.abundance.final$ID==i),]
  temp.sapro <- temp[which(temp$trophic.guild=="Saprotroph" & temp$Abundance > 0 & temp$confidence != "Possible"),]
  temp.symbio <- temp[which(temp$trophic.guild=="Symbiotroph" & temp$Abundance > 0 & temp$confidence != "Possible"),]
  temp.long <- temp[which(temp$distance.type=="long-distance" & temp$Abundance > 0),]
  temp.short1 <- temp[which(temp$distance.type=="short-distance_coarse" & temp$Abundance > 0),]
  temp.short2 <- temp[which(temp$distance.type=="short-distance_delicate" & temp$Abundance > 0),]
  temp.short3 <- temp[which(temp$distance.type=="contact" & temp$Abundance > 0),]
  temp.medium1 <- temp[which(temp$distance.type=="medium-distance_fringe" & temp$Abundance > 0),]
  temp.medium2 <- temp[which(temp$distance.type=="medium-distance_smooth" & temp$Abundance > 0),]
  temp.short.total <- rbind(temp.short1, temp.short2, temp.short3)
  exploration.guild.RA$saprotroph.richness[i] <- length(unique(temp.sapro$Species))
  exploration.guild.RA$symbiotroph.richness[i] <- length(unique(temp.symbio$Species))
  exploration.guild.RA$long.richness[i] <- length(unique(temp.long$Species))
  exploration.guild.RA$short.richness[i] <- length(unique(temp.short.total$Species))
  exploration.guild.RA$medium.fringe.richness[i] <- length(unique(temp.medium1$Species))
  exploration.guild.RA$medium.smooth.richness[i] <- length(unique(temp.medium2$Species))
}

#########################################################################
# DIVERSITY CALCULATIONS FOR TROPHIC GUILD AND EXPLORATION TYPE #
#########################################################################

# TROPHIC GUILD #

exploration.guild.RA$saprotroph.diversity <- c(1:96)
exploration.guild.RA$symbiotroph.diversity <- c(1:96)

for (i in 1:96){
  temp <- all.otu.abundance.final[which(all.otu.abundance.final$ID==i),]
  temp <- temp[which(temp$trophic.guild=="Saprotroph" & temp$Abundance > 0 & temp$confidence != "Possible"),]
  
  # remove NA from sample
  temp1 <- temp %>% drop_na(Species)
  
  # select all unique species
  species <- unique(temp1$Species)
  
  # create vector to put in proportion of total abundance for each species
  p <- c(1:length(unique(temp1$Species)))
  
  for (j in 1:length(p)){
    species.specific <- species[j]
    temp2 <- temp1[which(temp1$Species==species.specific),]
    p[j] <- sum(temp2$Abundance)/sum(temp$Abundance)
  }
  log.times <- c(1:length(unique(temp1$Species)))
  for (k in 1:length(p)){
    log.times[k] <- p[k]*ln(p[k])
  }
  exploration.guild.RA$saprotroph.diversity[i] <- -sum(log.times)
}

for (i in 1:96){
  temp <- all.otu.abundance.final[which(all.otu.abundance.final$ID==i),]
  temp <- temp[which(temp$trophic.guild=="Symbiotroph" & temp$Abundance > 0 & temp$confidence != "Possible"),]
  
  # remove NA from sample
  temp1 <- temp %>% drop_na(Species)
  
  # select all unique species
  species <- unique(temp1$Species)
  
  # create vector to put in proportion of total abundance for each species
  p <- c(1:length(unique(temp1$Species)))
  
  for (j in 1:length(p)){
    species.specific <- species[j]
    temp2 <- temp1[which(temp1$Species==species.specific),]
    p[j] <- sum(temp2$Abundance)/sum(temp$Abundance)
  }
  log.times <- c(1:length(unique(temp1$Species)))
  for (k in 1:length(p)){
    log.times[k] <- p[k]*ln(p[k])
  }
  exploration.guild.RA$symbiotroph.diversity[i] <- -sum(log.times)
}
  
# EXPLORATION TYPE #

exploration.guild.RA$long.diversity <- c(1:96)
exploration.guild.RA$medium.fringe.diversity <- c(1:96)
exploration.guild.RA$medium.smooth.diversity <- c(1:96)
exploration.guild.RA$short.diversity <- c(1:96)

for (i in 1:96){
  temp <- all.otu.abundance.final[which(all.otu.abundance.final$ID==i),]
  temp <- temp[which(temp$distance.type=="long-distance" & temp$Abundance > 0),]
  
  # remove NA from sample
  temp1 <- temp %>% drop_na(Species)
  
  # select all unique species
  species <- unique(temp1$Species)
  
  # create vector to put in proportion of total abundance for each species
  p <- c(1:length(unique(temp1$Species)))
  
  for (j in 1:length(p)){
    species.specific <- species[j]
    temp2 <- temp1[which(temp1$Species==species.specific),]
    p[j] <- sum(temp2$Abundance)/sum(temp$Abundance)
  }
  log.times <- c(1:length(unique(temp1$Species)))
  for (k in 1:length(p)){
    log.times[k] <- p[k]*ln(p[k])
  }
  exploration.guild.RA$long.diversity[i] <- -sum(log.times)
}

for (i in 1:96){
  temp <- all.otu.abundance.final[which(all.otu.abundance.final$ID==i),]
  temp.short1 <- temp[which(temp$distance.type=="short-distance_coarse" & temp$Abundance > 0),]
  temp.short2 <- temp[which(temp$distance.type=="short-distance_delicate" & temp$Abundance > 0),]
  temp.short3 <- temp[which(temp$distance.type=="contact" & temp$Abundance > 0),]
  temp <- rbind(temp.short1, temp.short2, temp.short3)
  
  # remove NA from sample
  temp1 <- temp %>% drop_na(Species)
  
  # select all unique species
  species <- unique(temp1$Species)
  
  # create vector to put in proportion of total abundance for each species
  p <- c(1:length(unique(temp1$Species)))
  
  for (j in 1:length(p)){
    species.specific <- species[j]
    temp2 <- temp1[which(temp1$Species==species.specific),]
    p[j] <- sum(temp2$Abundance)/sum(temp$Abundance)
  }
  log.times <- c(1:length(unique(temp1$Species)))
  for (k in 1:length(p)){
    log.times[k] <- p[k]*ln(p[k])
  }
  exploration.guild.RA$short.diversity[i] <- -sum(log.times)
}

for (i in 1:96){
  temp <- all.otu.abundance.final[which(all.otu.abundance.final$ID==i),]
  temp <- temp[which(temp$distance.type=="medium-distance_fringe" & temp$Abundance > 0),]
  
  # remove NA from sample
  temp1 <- temp %>% drop_na(Species)
  
  # select all unique species
  species <- unique(temp1$Species)
  
  # create vector to put in proportion of total abundance for each species
  p <- c(1:length(unique(temp1$Species)))
  
  for (j in 1:length(p)){
    species.specific <- species[j]
    temp2 <- temp1[which(temp1$Species==species.specific),]
    p[j] <- sum(temp2$Abundance)/sum(temp$Abundance)
  }
  log.times <- c(1:length(unique(temp1$Species)))
  for (k in 1:length(p)){
    log.times[k] <- p[k]*ln(p[k])
  }
  exploration.guild.RA$medium.fringe.diversity[i] <- -sum(log.times)
}

for (i in 1:96){
  temp <- all.otu.abundance.final[which(all.otu.abundance.final$ID==i),]
  temp <- temp[which(temp$distance.type=="medium-distance_smooth" & temp$Abundance > 0),]
  
  # remove NA from sample
  temp1 <- temp %>% drop_na(Species)
  
  # select all unique species
  species <- unique(temp1$Species)
  
  # create vector to put in proportion of total abundance for each species
  p <- c(1:length(unique(temp1$Species)))
  
  for (j in 1:length(p)){
    species.specific <- species[j]
    temp2 <- temp1[which(temp1$Species==species.specific),]
    p[j] <- sum(temp2$Abundance)/sum(temp$Abundance)
  }
  log.times <- c(1:length(unique(temp1$Species)))
  for (k in 1:length(p)){
    log.times[k] <- p[k]*ln(p[k])
  }
  exploration.guild.RA$medium.smooth.diversity[i] <- -sum(log.times)
}

#########################################
# incorporate C:N and other soil data #
#########################################

sample.data <- read.csv("soil.data.csv", header=TRUE)

# remove bog sample before using to input values 
bog.sample <- sample.data[which(sample.data$Transect == "Bog"),]
bog.RA <- exploration.guild.RA[which(exploration.guild.RA$Transect == "Bog"),]
sample.data <- sample.data[which(sample.data$Transect != "Bog"),]
exploration.guild.RA <- exploration.guild.RA[which(exploration.guild.RA$Transect != "Bog"),]

exploration.guild.RA$organic.soil.C.N <- c(1:95)
exploration.guild.RA$organic.soil.GWC <- c(1:95)
exploration.guild.RA$carcass.load <- c(1:95)
exploration.guild.RA$mineral.soil.C.N <- c(1:95)
exploration.guild.RA$GWC <- c(1:95)
exploration.guild.RA$carcass.deposition <- c(1:95)
exploration.guild.RA$decomposing.carcass <- c(1:95)
exploration.guild.RA$ph <- c(1:95)
exploration.guild.RA$slope <- c(1:95)
exploration.guild.RA$alder <- c(1:95)

# insert covariate data
for (i in 1:95){
  for (j in 1:94){
    if (exploration.guild.RA$Stream[i] == sample.data$Stream[j] & exploration.guild.RA$Distance[i] == sample.data$Distance[j] & 
        exploration.guild.RA$Side[i] == sample.data$Side[j] & exploration.guild.RA$Transect[i] == sample.data$Transect[j]){
      exploration.guild.RA$organic.soil.C.N[i] <- sample.data$organic.soil.C.N[j];
      exploration.guild.RA$carcass.load[i] <- sample.data$carcass.load[j];
      exploration.guild.RA$mineral.soil.C.N[i] <- sample.data$mineral.soil.C.N[j];
      exploration.guild.RA$organic.soil.GWC[i] <- sample.data$organic.soil.GWC[j];
      exploration.guild.RA$carcass.deposition[i] <- sample.data$carcass.deposition[j];
      exploration.guild.RA$decomposing.carcass[i] <- sample.data$decomposing.carcass[j];
      exploration.guild.RA$ph[i] <- sample.data$ph[j];
      exploration.guild.RA$slope[i] <- sample.data$slope[j]
      exploration.guild.RA$alder[i] <- sample.data$alder[j];
      }
  }
}

# put bog sample back in (although we don't really need to since we won't be using it)
bog.RA$C.N <- 39.03659652
bog.RA$GWC <- 0.66503545
bog.RA$carcass.deposition <- "no"
bog.RA$decomposing.carcass <- "yes"
bog.RA$ph <- "none"
exploration.guild.RA <- rbind(exploration.guild.RA, bog.RA)
sample.data <- rbind(sample.data, bog.sample)

write.csv(exploration.guild.RA, "exploration.guild.RA.csv")

